arXiv:1505.02962vl [math.NA] 12 May 2015 


Discrete line integral method for the Lorentz force system 


Haochen Li°, Yushun Wang a ’* 

a Jiangsu Key Laboratory for NSLSCS, School of Mathematical Sciences, Nanjing Normal University, 

Jiangsu, 210023, People’s Republic of China 


Abstract 

In this paper, we apply the Boole discrete line integral to solve the Lorentz force system 
which is written as a non-canonical Hamiltonian system. The method is exactly energy- 
conserving for polynomial Hamiltonians of degree v < 4. In any other case, the energy 
can also be conserved approximatively. With comparison to well-used Boris method, 
numerical experiments are presented to demonstrate the energy-preserving property of 
the method. 
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1. Introduction 


Geometric numerical integration methods have come to the fore, partly as an alter¬ 
native to traditional methods such as Runge-Kutta methods. A numerical method is 
called geometric if it preserves one or more physical/geometric properties of the system 
exactly (i.e. up to round-off error) 0,3 Si- Examples of such geometric properties 
that can be preserved are first integrals, symplectic structures, symmetries and reversing 
symmetries, phase space volumes, Lyapunov functions, foliations, etc. Geometric meth¬ 
ods have applications in many areas of physics, including celestial mechanics, particle 
accelerators, molecular dynamics, fluid dynamics, pattern formation, plasma physics, 
reaction-diffusion equations, and meteorology 0, [(3, 00, $M- 


Many important phenomena in plasmas can be understood and analyzed in terms of 
the single-particle motion which satisfies the Lorentz force equations [ll|. The motion 
of charged particles in single particle model is governed by the Newton equation under 
the Lorentz force exerted by a given electromagnetic field. Hamiltonian formulation is 
available for the single particle model and other magnetized plasma models in special co¬ 
ordinates 12,113]. In long-term simulating the motion of charged particles, non-geometric 


methods such as the standard 4th order Runge-Kutta method may rise to a complete 
wrong solution orbit, since the numerical errors of each time step will add up coherently 
and become significantly large over many time steps. In contrast, geometric methods 
show a good long-term accuracy, for instance, the volume-preserving algorithms (ljj |. the 


Boris method which can also conserve phase space volume 14, [15|, [16], [IT], the varia¬ 
tional symplectic method 19, 2c], 21] and so on. It is well-know that the most noticeable 
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23]. The averaged vector field method 


structure of a Hamiltonian system is the Hamiltonian function itself which is usually the 
energy of the system. The paper is devoted to constructing an energy-preserving method 
for the system. 

The conservation of the energy function is one of the most relevant features char¬ 
acterizing a Hamiltonian system. Methods that exactly preserve energy have been 
considered since several decades. Many energy-preserving methods have been pro¬ 
posed [ 22 I, 23, 0, 0, 26]. The discrete gradient method is among the most popular 
methods for designing integral preserving schemes for ordinary differential equations, 
which was perhaps first discussed by O. Gonzalez [ 22 }. T. Matsuo proposed discrete 
variational method for nonlinear wave equation 
which is a B-series method has been proposed [2 


25| . More recently, L. Brugnano and 
F. Iavernaro proposed the discrete line integral (DLI) methods (0 . 28 ] and the Hamil¬ 
tonian boundary value methods 29, 3o|. The key idea of the DLI methods is to exploit 
the relation between the method itself and the discrete line integral, i.e., the discrete 
counterpart of the line integral in conservative vector fields. This tool yields exact con¬ 
servation for polynomial Hamiltonians of arbitrarily high-degree. Different qudarature 
formulas yield different DLI methods. Especially, if we choose the Boole’s rule which is 
the Newton-Cotes formula of degree 4, the so-called Boole discrete line integral (BDLI) 
method which is exactly energy-conserving for a polynomial Hamiltonian of degree v < 4 
is obtained. 

In this paper, the Lorentz force system is written as a non-canonical Hamiltonian 
form. We apply the BDLI method for the Hamiltonian system, and a new energy¬ 
preserving method is obatined. The new method is symmetric and can preserve the 
Hamiltonian up to round-off error. 

RefBoris 

The paper is organized as follows: In Sect. [2l the dynamics of charged particles in the 
electromagnetic field is shown and it is written as a non-canonical Hamiltonian system. 
In Sect. El we use the BDLI method to solve the Hamiltonian system. Based on the 
Boole’s rule, a new method for Lorentz force system is obtained. With comparison to 
well-known Boris method 14], numerical experiments are presented in Sect. [4] to confirm 
the theoretical results. We finish the paper with conclusions in Sect. [5j 


2. Hamiltonian form of the Lorentz force system 


In this section, we review the Hamiltonian form of the Lorentz force system [12j . 1131 . 


18] - For a charged particle in the electromagnetic field, its dynamics is governed by the 


Newton-Lorentz equation 


mx = g(E + x x B), x G M 3 , (2.1) 

where x is the position of the charged particle, m is the mass, and q denotes the electric 
charge. For convenience, here we assume that B and E are static, thus B = V x A and 
E = — V(/3 with A and the potentials. 

Let the conjugate momentum be p = mx + gA(x), then the system (12.ip is Hamil¬ 
tonian with 

#(x, P) = tt-(p “ ? A (x)) ■ (P - ?A(x)) + g<p(x). (2.2) 
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(x, v),x = x, v = p/m — 

(2.3) 

(2.4) 


Recasting (12 .1 1) with transformation G : (x, p) — ; 
gA(x)/m, we can obtain 

x = v, 

V = — (E(x) + V x B(x)). 
m 

Denote z = [x T , v I ] r . (12.31b(12.4|) can be written as a non-canonical Hamiltonian system 

z = /(z) = K(z)X7H(z), (2.5) 

where H(z ) = mv • v/2 + qip(x), and 

0 


K( z) = 

is a skew-symmetric matrix with 


±1 

17 % 

-±I -%B(x) 

m m z k / 


0 B 3 (x) B 2 (x) 

B(x) = | —B 3 (x) 0 £?i(x) 

D 2 (x) -B i(x) 0 


defined by B(x) = [Bi(x), R 2 (x), -B 3 (x)F 


3. Boole discrete line integral method 


It is well know that the flow of the system (12.5h preserves the energy which is usually 
the Hamiltonian H( z) exactly. In this section, we derive a new en ergy -preserving scheme 
for the system (12.51) by using the BDLI method proposed in 271, [28]. Starting from the 
initial condition z 0 we want to produce a new approximation at t = h, say z 1; such that 
the Hamiltonian is conserved. By considering the simplest possible path joining z 0 and 
zi, i.e., the segment 

<j(ch) = czi + (1 — c)z 0 , c G [0,1], 

we obtain that 


kff( zi) - H(z „)) = i(tfMft)) - H(a( 0))) 


h 


h 

1 


(3.1) 

(3.2) 


— — I \/H(a(t)) T a'(t)dt 


h 


= / VH(cr(ch)) o\ch)dc 

Jo 

1 f 1 

= h J V H(czi + (1 - c)z 0 ) T (zi - z 0 )dc 
= [ [ X7H(czi + (1 - c)z 0 )dc] T Zl - — = 0, 


h 


provided that 


zi — z 0 

h 


= K (z) f VH(cz 1 + (1 — c)z 0 )dc, 
Jo 


(3.3) 
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where we choose z = Z1 + Z0 in this paper. In fact, due to the fact that K ( Zl ^ Z(1 ) is skew 
symmetric, we have 

[J V/f(cz! + (1 - c)z 0 )dc] T Zl — Z ° (3.4) 

= [ f VH(czi + (1 - c)z 0 )dc] T K( — ^ Z ° )[ f WH(czi + (1 - c)z 0 )dc] = 0. 

Jo * Jo 

In this paper, we use the Boole’s rule to calculate the integrand at the right-hand 
side in (j3.3(1 numerically. Then we obtain the BDLI method for (12.51) : 

Zl =z„ + Lk(^L±^ i)(7Vff(z 0 ) + 32Vff (iR±£i) (3.5) 

+ 12 Vff( Z °^ Zl ) + 32Vg( Z °^ 3zi ) + 7Vff(z,)). 

The method (j3.5[) is exactly energy-conserving if H is a polynomial Hamiltonian of degree 
v < 4. It is obvious that the method is symmetric and, therefore, of order 2. By the 
mono-implicit character of the method, the computational cost is less than high order 
implicit Runge-Kutta methods which need the values of every stages. 

Remark 3.1. Actually, we can also use other quadrature formulas to calculate the in¬ 
tegrand numerically, for instance, the trapezoid rule, the Simpson’s rule and so on. We 
observe that for any non-polynomial Hamiltonian, energy conservation can be practically 
obtained, by choosing a suitable quadrature formula such that the errors of the Hamil¬ 
tonian are up to the round-off error 28]. In fact, from the numerical experiments we 
know that the non-polynomial Hamiltonian can be preserved up to the round-off error of 
1(T 15 if the Boole’s rule is used (see Fig. 2(b)). 

4. Numerical experiments 

In this section, we numerically test the BDLI method (13.51) . The fixed point iterative 
method is used. 

4.1. 2D dynamics in a static electromagnetic field 

Firstly, we consider a non-polynomial Hamiltonian case, i.e., the 2D dynamics of the 
charged particle in a static, non-uniform electromagnetic field 

1CT 2 

B = V x A = Re z , E = — = -^-(xe x + ye y ), (4.1) 

where the potentials are chosen to be A = <p = hLi j n cylindrical coordinates 

(R,£,z) with R = \Jx 1 + y 2 . In this example, the physical quantities are normalized by 
the system size a, the characteristic magnetic field Bo. It is known that the energy 

1 1(R 2 

R=-v.v+— (4.2) 

is a constant of motion. As the given electromagnetic held (14.Tj) changes slowly with 

v 2 

respect to the spatial period of the motion, the magnetic moment /1 = ^ is an adiabatic 
invariant, where Vi is the component of v perpendicular to B. 
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Figure 1: The Boris method is applied to the simple 2D dynamics with step h = 7r/10. (a) The orbit 
in the 100-st turn; (b) Errors of the angular momentum p the magnetic moment p and the energy H 
for t £ [0, 5 x 10 4 h]. 



(a) x (b) time 

Figure 2: The BDLI method is applied to the simple 2D dynamics with step h = 7r/10. (a) The orbit 
in the 100-st turn; (b) Errors of the angular momentum p £, the magnetic moment p and the energy H 
for t G [0, 5 x 10 i h\. 
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Starting from the initial position Xo = [0,0.1,0] r with the initial velocity Vo = 
[0.1, 0.01, 0] T , the analytic orbit of the charged particle is a spiraling circle with constant 
radius. The large circle corresponds to the V • B drift and the ExB drift of the guiding 
center, and the small circle is the fast-scale gyro-motion. 

We first apply the Boris method and the numerical results are shown in Fig. 1. The 
step size is h = 7r/10 which is the 1/20 of the characteristic gyro-period. Fig. 1(a) shows 
the numerical orbits of the Boris method in the 100-st turn. In Fig. 1(b), we illustrate 
the errors of the angular momentum p g, the magnetic moment p and the energy H for 
t e [0, 5 x 10 4 /i]. It is observed that the errors are bounded over a long integration time. 

Next we text the BDLI method (13.51) with the same step and the numerical results 
are shown in Fig. 2. Fig. 2(a) shows the numerical orbits of the BDLI method in the 
100-st turn. Fig. 2(b) shows errors of the invariants for t G [0,5 x 10 4 /r]. It can be 
observed that the errors of the angular momentum p^ and the magnetic moment p are 
also bounded over a long integration time. What is more, the BDLI method can preserve 
the non-polynomial Hamiltonian (14.2[) up to the round-off error 


31 


4-2. 2D dynamics in an axisymmetric tokamak geometry 

In this subsection, we consider a polynomial Hamiltonian case, i.e., the motion of a 
charged particle in a 2-dimensional axisymmetric tokamak geometry without inductive 
electric field. The magnetic field in the toroidal coordinates (r, #,£) is expressed as 


B 


B 0 r 

qR 


e 9 + 


B(\ Ho 

R 




(4.3) 


where B 0 — 1, R 0 — 1, and q = 2 are constant with their usual meanings. The corre¬ 
sponding vector potential A is chosen to be 


A 


z (1 -R) 2 + z 2 In R 

e R H-tt; - e £ H-— e z- 


2 R 


4 R 


(4.4) 


In this example, the physical quantities are normalized by the system size a, characteristic 
magnetic field B 0 . In the absence of electric field, the energy of the system H — |v ■ v is 
a polynomial Hamiltonian. With the conserved quantities, the solution orbit projected 
on (R, z ) space forms a closed orbit. 

To apply the discrete line integral method, we transform B (14.3[) into the Cartesian 
coordinates (x, y, z ) which is 


2 y + xz 2 x — yz R — 1 

2 R 2 2 R 2 y 2 R y 


(4.5) 


Starting with the initial position x 0 = [1.05, 0, 0] T and the initial velocity v 0 = [0,4.816 x 
10 -4 , 2.059 x 10~ 3 ] T , the orbit projected on (R,z) space is a banana orbit, and it will 
turn into a transit orbit when the initial velocity is changed to Vo = [0,2 x 4.816 x 
ICC 4 , 2.059 x 10' 3 ] T . Setting the step size h = 7r/10 which is the 1/20 of the gyro-period, 
we adopt the BDLI method. The simulation over 5 x 10 4 steps is shown in Fig. 3. Fig. 
3(a) shows the banana orbit solution and Fig. 3(a) shows the transit orbit solution. It 
is observed that the BDLI method gives correct orbits over a very long integration time. 
In the case without an electric held, the BDLI method can also preserve the energy H 
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Figure 3: Numerical solution of the BDLI method with step h = 7r/10 for t G [0, 5 x 10 5 h]. (a) Banana 
orbit; (b) Transit orbit; (c) Energy preservation. 


exactly. In Fig. 3(c) we display the errors of the energy. It can be observed that the 
errors of the energy are up to the round-off error of ICC 18 . Due to the iterative method 
which is not energy-preserving, a drift is present in the numerical Hamiltonian. The 
above experiments show the superiority of the BDLI method for solving the Lorentz 
force system with a static electromagnetic held. 

5. Conclusions 

In this paper, we propose the BDLI method for the Lorentz force system which is 
written as a non-canonical Hamiltonian system in the coordinates (x, v). The numerical 
results show that, for polynomial Hamiltonian case, the method can conserve the energy 
exactly, and for non-polynomial Hamiltonian case, the method can also conserve the 
energy up to round-off error. Due to the energy-preserving property, the BDLI method 
has good performance over long simulation time. 
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